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ABSTRACT 

We test the accuracy of different models of the attenuation of light due to resonant scattering 
by intergalactic neutral hydrogen by comparing their predictions of the evolution of the mean 
cosmic flux decrement. Da, to measurements of this quantity based on observations. To this 
end, we use data available in the literature and our own measurements of the cosmic flux 
decrement for 25 quasars in the redshift range 2.71 < Zem < 5.41 taken from the SDSS Data 
Release 5. In order to perform the measurements of Da, we fit a power-law to the continuum 
redward of the Lyff emission line, and extrapolate this fit to region blueward of it, where the 
flux is severely affected by absorption due to intervening Hi absorbers. 

We compute, using numerical simulations, the redshift evolution of Da accounting for 
the presence of Lya Forest absorbers and Lyman limit systems randomly distributed along the 
line-of-sight, and compute its intrinsic scatter at the 1-, 2-, and 3 cr level due to fluctuations 
in the absorber properties (column density, Doppler parameter, redshift) along different lines- 
of-sight. The numerical simulations consist of Monte Carlo realizations of distributions of the 
absorber properties constrained from observations. 

The results from the models considered here confirm our theoretical expectation that the 
distribution of D a at any given redshift be well described by a lognormal distribution function. 
This implies that the effective optical depth, usually defined as the negative logarithm of the 
average flux, 1 - Da, is very accurately Gaussian distributed, in contrast to previous studies. 
This result is independent to the form of the input distribution functions, and rather insensitive 
to the presence of high-column density absorbers, such as the Lyman limit systems. 

By comparing our and previous measurements of Da to the outcomes of our simulations, 
we find an excellent agreement between the observations and the evolution of the mean D f^?& 
predicted by one of the models considered in this work. The observed scatter in Da at each 
redshift, however, cannot be recovered from our simulations. Even though there is evidence 
for the fact that the lack of agreement between models and observations comes from the 
combination of heterogeneous measurement sets obtained by different methods, the failure of 
the models to accurately account for the absorption by intergalactic Lya absorbing systems 
and its variation along different lines-of-sight cannot be completely ruled out. 

Key words: methods: numerical, intergalactic medium, quasars: absorption lines 



1 INTRODUCTION 

Since the introduction of the Gunn-Peterson (GP) test by 
iGunn & PetersonI l ll965h . a detailed knowledge about the physical 
state of the intergalactic medium (IGM) has been gained from the 
study of the absorption features identified in the spectra of quasi- 
stellar objects (QSOs) at restframe wavelengths A < 121.5 nm, 
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which are now known to be mainly due to resonant scattering by 
intergalactic neutral hydroge n randomly di stributed along the line- 
of-sight, as first proposed by iLvnd j jl97lh . For instance, the null 
result in the search for a GP trough has bee n used to rule out 
the existence of a hot intercloud m edium (ICM) l lSteidel & Sargenj 
ll987allGiallongo et al.|[l992l.[l994t), which was thought to confine 
by pr essure the Lya clouds Isargent et al.ll98 0. Ostriker & Ikeuchil 
Il983h . As a result of detailed analyses of the line statistics of 
the abso rbing material, a we alth of information on its clustering 
(see e.g. lOstriker et al.lll98^ , in particular the existence of voids 
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Pierre et al.|[l988l , lOuncan et al1ll989l. iDobrzvcki & Bechtold 



I199I ), and the evolution of its number de nsities, column d ensities, 
and Doppler parameter with redshift (e. e. iKim et al ] l 19971) has ac- 
cumulated over the past years. These results, in combination with 
the use of state-of-the-art numerical simulations of structure forma- 
tion based o n the currently accep ted paradigm of the concordance 
cosmology dSpringel et alj 12003) . show that the features seen in 
absorption against bright background sources arise when the line- 
of-sight intersects the structures that naturally emerge and evolve 
with time under the influence of gravitational attraction. Differ- 
ent types of structures such as the filaments present in the cosmic 
web, galactic haloes, and even the discs of primeval galaxies, give 
rise to distinct absorption features attributed to entities known his- 
torically as Lyo- Forest clouds, Lyman l imit systems (LLSs), and 
damp ed Lya absorbers (DLAs) (see e.g. lRauchlll998llWolfe et al] 
l200i for excellent reviews, respectively). Nevertheless, the rela- 
tion between the observed absorption features and the objects caus- 
ing them, in particular the correlation between observed damped 
absorption lines, metal lines (e.g. Mg n, O vi) and galaxies-the so- 
call ed Absorber-galaxy Connection -is still a matter of debate (see 
e.g. IWilliams, Shu & Menard [200I Part 1). As a consequence of 
numerous efforts over many years, we now have a better under- 
standing of the origins of the different absorption features observed 
in QSO spectra. In particular, the notion of discrete, intervening Hi 
absorbing systems randomly distributed along the line-of-sight has 
been embedded into the more general picture of an evolving con- 
tinuous intergalactic medium with a Hi density field that varies in 
space and time, with its evolution driven mainly by the Hubble ex- 
pansion, the radiation field of ionising UV sources, and the collapse 
of structures due to gravity. 



1.1 Methods and Input Distributions: A Brief Review 

Over decades many people have been working hard towards in- 
ferring the physical properties of the intergalactic medi um such as 
its ch emical content, density, temperature, etc. (see e.g. I Kim et al] 
l2002i) by measuring the type of transition, strength, number den- 
sity, and profiles of absorption lines imprinted in the spectra of 
QSOs and, more rece ntly, of Gamma-Ray bursts (GRBs) (e.g. 
iLamb & Reichartl2000l) . There has also been a great effort to quan- 
tify the effect of the absorption due to intergalactic neutral hy- 
drogen on the photometric properties of background sources. As 
a matter of fact, several models have been developed in order to 
account for this so-called intergal actic attenu ation, with different 
approaches and purposes. M0ller & JakobsenI (1990) used Monte 
Carlo simulations to estimate the amount of absorption at wave- 
lengths shorter than the redshifted Hen A 30.4 nm line, in order to 
test the feasibility of the equivalent of the Gunn-Peterson test for 
intergalactic helium. They found that the absorption as a function 
of wavelength, averaged over many lines-of-sight, should display 
together with a characteristic stair-case profile due to the cumula- 
tive absorption at the Hi resonant wavelengths, an additional char- 
acteristic valley-shaped feature (the "Lyman-valley") due to the cu- 
mulative effect of the photoionisation of Hi by photons with ener- 
gies Ey > h c /Ai^, where = 91.2 nm. Later on, in a seminal 
paper iMadad j 19951) developed an analytical method to quantify 
the opacity due to intergalactic Hi as a function of redshift, and 
its effect on the colors of high-redshift galaxies. The underlying 
assumption of this model is that the observed flux of a source at 
redshift z is given by the product of the intrinsic flux and a trans- 
mission factor that accounts for the mean absorption as a func- 
tion of wavelength given in the form exp(-Te(f) = (exp(-T)) (see 



lMadadll995l , their equation 3), where the brackets denote the aver- 
age over an ensemble of random lines-of-sight. The most common 
application of this model consists in correcting the flux of a syn- 
thetic (galaxy-, QSO-) spectrum for intergalactic absorption. This 
correction is of particular importance at high redshift, where in- 
tergalactic Hi severely absorbs the light of a background object at 
restframe wavelengths shorter than 121. 6nm, leading t o a substan- 
tial reddening of its colour (see e.g. iBicker e t al."2004'). As the nu- 
merous references in the literature attest, the Madau model has be- 
come tiiejnost_wddel;^jised_attem However, in a later 
work lBershadv. Charlton & GeofTrovl h999l . from now on BCG99) 
argued that it is not possible to estimate the mean change in the 
magnitude of a source at a given redshift due to absorption by in- 
tergalactic Hi along the line-of-sight by multiplying the mean trans- 
mission curve of IMadad ' s model with the spectrum of the source 
and integrating over the corresponding passband, mainly because of 
the existence of color terms. They suggested that the correct way 
of accounting for the mean effect of Hi absorption on the spec- 
trum of a background source and on its photometric properties is 
to model first the absorption along many random lines-of-sight, 
compute the desired photometric quantities for each one of them, 
and then compute the average over the ensemble of lines-of-sight. 
In other words, they argue that the processes of averaging over 
many random lines-of-sight and measuring photometric quantities 
are non-commutative. Indeed, they showed using a Monte Carlo 
technique that the average magnitudes computed follo wing the ir 
approach substantially differ from those computed using iMadad 's 
model, even when using the same input distributions for the number 
of absorbers, their col umn den sities and Doppler parameters. The 
approach proposed by lBCG99l efl^ectivelv mimics the measurement 
process that would take place if one would determine e.g. the mean 
observed brightness of a collection of galaxies with different ab- 
sorber populations along their particular lines-of-sight, but other- 
wise identical in their intrinsic properties (spectrum, morphology, 
etc.), and is hence physically meaningful. It turns out that features 
such as the characteristic stair-case profile and the Lyman-valley 
cannot possibly be observed in a single spectrum, since they arise 
only by averaging over sufficient numbers of lines-of-sight, a pro- 
cess that has no physical me aning. 

In a more recent paper, iMeiksid j2006l) developed a method 
to compute the opacity due to intergalactic Hi by using hydrody- 
namical simulations of structure formation in the f ramework of the 
conco rdance ACDM cosmology performed by IMeiksin & White! 
( l2004h . Applying their model to compute broad-band magnitudes 
for differ ent types of obj ects (e.g. starburst galaxies, QSOs of Type 
I and II), Meiksin !200&) reports differences of 0.5 - 1.0 mag with 
respect to Madau ( 1995)'s mo del. Despite th e different results ob- 
tained, this model is similar to lMadad jl995l) 's model in the sense 
that it implicitly assumes that the mean opacity of the IGM along 
a random line-of-sight due to the presence of Hi can be accounted 
for by multiplying a given input spectrum with a mean attenuation 
curve of the form exp( -Teff) and i ntegrating over the corresponding 
filter function (see M eiksidl2006l their equation 8). 

Following iBCG99l we state that 



1: 



/,r<exp(-T))-r(i)d^# 



f,-exp(-T)-T(A)dA) ,(I) 



where is the intrinsic flux, T(A) is the filter transmission func- 
tion, and the brackets denote the average over all lines-of-sight. 
We consider that the operation implied by the right-hand side of 
this expression is the correct way of estimating mean magnitudes 
of background objects including the effect of the absorption due 
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to intergalactic Hi. This approach is of course not restricted to the 
computation of mean magnitudes and colors, and can be applied to 
the estimate of the mean of any photometric quantity. Furthermore, 
it is also possible to determine not only the mean, but in principle 
any desired confidence interval around the mean, e.g. ±cr range, 
via the computation of quantiles (see Section|5]l- 

It should be clear that not only the method, but also the input 
physics is an (even more) crucial ingredient of a particular model 
that a ccounts for the intergalactic attenuation, as already shown by 
IbCG 99. It is, however, not trivial to test whether using a particular 
method and a set of input distributions accurately describes the 
observed effect of the absorption by intergalactic Hi on the spectra 
of background sources. For example, the evolutionary synthesis 
models of Bicker etal. (2004) that include the correction for 
intergalactic absorption based on Madau's model mat ch quite well 
the ob servations of galaxies in the Hubble Deep Field l lBicker et al.l 
|2004 their Figure 12), since the magnitude differences reported 
byjBCG99 (their Figure 7) for one of their models with respect 
to|Madau's model are in this case of the order of the scatter of 
the observations around the predicted colors. In other words, 
even though iMadaui 's and iBCG99i 's approaches are fundamentally 
different, it is difficult to test the accuracy of their predictions on 
the basis of a comparison to e.g. observed galaxy colors. A quantity 
that is more sensitive to the absorption due to intergalactic Hi is the 
mean cosmic flux decrement Da (cf. Section |2}- The reason for 
this is that the restframe wavelength range over which this quantity 
is measured is typically 10 nm wide, and is hence narrower than 
typical broadband filters. We thus consider as a primer test that 
any model that accounts for the absorption due to intergalactic 
Hi should reproduce first of all the observations of this quantity. 
Hence, through the comparison of their respective predictions to 
measurements of Da it should be possible to discriminate between 
different models which may or may not be appropriate to account 
for the effect of the intergalactic attenuation on the spectra of 
background sources. We select for this comparison two of several 
models introduced bv lBCG99l to compute the magnitude changes 
of high-redshift galaxies due to intergalactic absorption, and 
which ar e ba s ed, re spectively, on the input distrib ution functions 
of .Kim et^ h997h . and the input physics of the iMadaJ j 1 9951) 
approach. The latter is chosen since, as already mentioned, it is the 
model most widely used in the literature; the former constitutes, 
to the best of our knowledge, the most complete set of input 
distribution functions for the evolution of the hya ab sorbers to 
date, later expanded but not significantly changed by I Kim et al] 
( I2OO2I) . 

To sum up: The main goal of this work is to model and analyse 
the redshift evolution of Da for different evolution scenarios of 
the intergalactic neutral hydrogen, conveniently parametrized by 
input distribution functions of the form of equation ^ Section^. 
Hence, the models we consider here differ only by the set of input 
distribution functions used, but they all equal in method, i.e. all 
of them take advantage of the Monte Carlo technique. We judge 
the goodness of a particular model by its power to reproduce the 
observations of Da in a wide redshift range, and analyse to which 
extent the theoretical expectations about the properties of Da are 
recovered from the measurements. 

This work is organised as follows: In Section |2] we briefly re- 
call the concept of the cosmic flux decrement and discuss some 
issues related to its measurement. In Section[3]we present two dif- 
ferent types of models for the intergalactic attenuation, which we 
use to model the redshift evolution of the cosmic flux decrement. 



In Section 151 we describe our measurements of this quantity for a 
sample of SDSS QSO spectra. Finally, we compare these and pre- 
vious measurements to the outcomes of each model, and discuss 
the results of this comparison as well as some other implications of 
the models for the evolution of Da in Section|5] 



2 THE COSMIC FLUX DECREMENT REVISITED 

Before high-resolution (i.e. AA < I nm), high S/N observation be- 
came feasible, the basic spectroscopic technique used to analyse the 
effect of the absorption due to intergalactic neutral hydrogen on the 
spectra of background sources was to measure the mean depression 
of the observed flux relative to the unabsorbed flux -or emission 
continuum-, a quantity which became to be known as cos mic flux 
decrement. This quantity, first introduced by Oke & Korvcanskvl 
( Il982h . can be defined as a function of redshift by 

where and are the continuum and the observed fluxes, re- 
spectively, and AA = (I + z) ■ (A2 - A i). Formally, the integral is 
computed in the restframe wavelength range [102.5, 121.6] nm, i.e. 
between the Ly^S and the Lyor emission lines. However, the actual 
estimate of Da is usually performed between the restframe wave- 
lengths /I = 105 nm and A = 111 nm -or in an even narrower wave- 
length interval- in order to avoid contamination by the emission 
wings of the Ly/J + Ovi and Lya lines, respectively. 

Since Da is extremely sensitive to f^, as can be easily seen 
from the definition ([2j, an accurate measurement of this quantity 
demands a reliable estimate of the underlying continuum. Unfor- 
tunately, there is no consensus of what the best method to esti- 
mate the continuum may be. A popular choice, mainly because 
of the presence of emission lines redward of the Lyor forest re- 
gion, consists in fitti ng a l ocal continuum, most commonly using 
cubic splines (see e.^ . lLu et al. 1996) or b-spline functions (see e.g. 
iKirkman et alj|2003l iTvtler et al . 2004a b), searching for regions 
apparently free of absorption blueward of the Lyo- emission line. 
Other authors prefer to fit a continuum in the region redward of the 
red wing of the hya emissi on line, and ext r apolate it to the region 
blueward of it (^g ^j del & SargentI Il987bl . ISchneider et al] 1 19891 
ICristiani et al.ll 19931) . A widely adopted form for the fitted contin- 
uum in this case is a power-law with spectral index a y which for 
wavelengths A > 121.6 nm takes on values in the ran ge [0.28, 0.99] 
(Steidel & Sargentlll987bl . lVanden Berk et aI.ll200lL and references 
therein). The latter method may tend to place the intrinsic contin- 
uum level higher than it actually is, thus overestimating the m ea- 
sured values of Da (cf. Sec. 14.11 see also lTvtler et al.ll2004 a'. and 
references therei n); for the former method th e opposite is true, in 
general (see e.g. iFaucher-Giguere et al. I I2OO7I) . For either method, 
there is an uncertainty in the estimate of the continuum, and this 
is the main drawback of the mean flux depression as a technique 
to estimate the mean absorption due to neutral hydrogen present 
in the int ergalactic m edium (IGM). In an attempt to overcome this 
problem. lZuo & ^(19931) used the idea that the cosmic flux decre- 
ment effectively measures the total equivalent width of all Lyo- ab- 
sorption lines in the chosen wavelength range -if corrected for the 
contribution of metal absorption lines- to measure this quantity by 
simply adding up the equivalent widths of lines identified as hya 
absorption lines. Of course, the reliability of this technique highly 
depends on the correct identification of lines, and misidentification 
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may lead to larger uncertainties than those associated to the contin- 
uum estimate. 

It turns out that reliable measurements of £)a are very use- 
ful to constrain estimates of fundamental cosmological parame- 
ters suc h as the mean bary on density fl^, the UV background in- 
tensity ( iRauch et al.lll99'^^ . the normalization of the power spec- 
trum CTg the vacuum-energy density Qa. and the Hubble param- 
eter Ho ( lTvtleretalJl2004 j). It is hence of great interest to con- 
trast observations to theoretical models of the evolution of Z) a, in 
which the bias due to the uncertainty in the estimate of the contin- 
uum is a bsent. This has been previously done by dilferent workers 
(see e.g. iGTallongo, Gratton & Treveseli 199(1 ICristiani et al.lll993l 
lMadatj|l995h T usually obtaining a good agreement with observa- 
tions. However, there is still a scatter in the observations of this 
quantity for which it has not been accounted yet in any modelling 
so far. It is possible that this scatter may well be due to the dif- 
ference in the methods used by each group to measure Da and to 
the different redshift ranges probed. Hence, by taking advantage of 
the Monte Carlo technique, we assess to which extent the obser-ved 
scatter can be ascribed to the intrinsic scatter in Da due to fluctu- 
ations in the properties (number density, column density, Doppler 
parameter) of the Hi absorbers along different lines-of-sight (see 
Section[33]l. 

From the theoretical point of view, assuming that the restframe 
equivalent width of the absorbers does not evolve with redshift, 
and that the number density of the absorbing systems evolves like 
oc (1 + z)'', it is expected that Z)a should evolve with z like 

£»a(z)'x(1+z)'+'' (3) 

where the extra factor (1 -I- z) comes from the scaling of the equiva- 
lent width, as pointed out by Jenkins &. Ostriker ( 1991). Indeed, it 
has been found empirically that the redshift evolution of £)a can be 
described by a power law Da (z) = A • (I -t-z)'^'' with A = 6.2-10"' 
and y = 1.75 jKirkman et aflllOO^ . even though other functional 
forms, e.g. an exponential of the form Da (z) = D^ - e"''*'' with 
D° = O.OI and a = 0.75 (Zhan g et al. .19971) a lso m atch well 
the observations. More recently, Kirkman et al. I ( l2007l) showed 
that the observed evolution of Da with redshift in the range 
< z < 3.2 is well described by a broken power-law, even though 
the significance of the fit is low. In any case, expressions like 
these are only valid up to a given redshift, since they diverge for 
z — > CO, whilst Da converges asymptotically to 1 in this limit 
or, more precisely, when z approaches the redshift Zreion at which 
reionisation sets on. As will be shown later, the redshift evolution 
of Da predicted by our simulations does satisfy this asymptotic 
behaviour (cf. Section|5j. Furthermore, even if the power-law form 
for the evolution of Da holds, the index y in equation ([3]l should 
be replaced by y, where the latter index accounts for the averaged 
evolution of absorbers of different column densities, which evolve 
all at different rates. Conversely, estimates of a single y from Da 
measurem ents assuming a pow er-law of the form of equation lO, 
as done bv lO'Brien et alj ( Il988l) . may give a hint on the population 
of absorbers dominating the behaviour of Da, comparing the 
estimated y with the power-law index of the different populations. 
For instance, according to one of the models considered here, 
the redshift evolution of Da will be shown to be dominated by 
absorbers with column densities Nhi < 10'^ cm"-. 

A compilation of Da measurements, accumulated in the liter- 
ature over the past twenty years approximately, and which includes 
our own measurements that extend the redshift range to z™, = 5.41, 



Table 1. Typ e s of a bsorbers and their corresponding parameters adopted 
from iMadaul j 19951 their equation 10). Note, however, that they quote 
Nhi = 2.0 - 10 as lowest column density, while we use Nhi = 1.0-10 
in order for the adopted normalisation to be consistent. 



Nhi [cm -2] 




No 


y 


P 


10'2 - 1.59- 


10" 


2.40 - lO' 


2.46 


1.50 


1.59- IQi' - 


1020 


1.90- lO** 


0.68 


1.50 



is shown in Figure [6] (cf. Section |4] and cited references for details 
on the measurements in this figure). 



3 MODELLING THE INTERGALACTIC ATTENUATION 

Since the observation of individual sources (galaxies, QSOs, 
GRBs) necessarily implies observations along different lines-of- 
sight, it is expected that the stochastic nature of the distribution of 
the Lya absorbers, especially of those with the highest column den- 
sities, causes a scatter in the observed absorption, even for sources 
with identical intrinsic spectra. Hence, depending on the absorp- 
tion along a particular line-of-sight, one would expect different ob- 
served values for each measurement of any photometric quantity, 
for example, the cosmic flux decrement Da. Performing enough 
measurements of such a quantity for sources with -ideally- iden- 
tical SEDs at a fixed redshift, one could in principle estimate its 
mean and its scatter due to stochastic effects in the absorption by 
neutral hydrogen in the IGM. 

The numerical realisation of this thought experiment is 
best achieved throu gh M onte Carlo simulations. Following e.g. 
'M0ller & Jakobsen' (1990'). lGiallongo et al.l l ll990l) . ICristiani et al.l 
( 1993), and BCG99, we generate thousands (4 • 10') of lines-of- 
sight each with a random population of Hi absorbers, and compute 
the absorption along each of them for a given input spectrum at 
a fixed redshift. The population of each line-of-sight consists of a 
random number Woj, of absorbing systems, each of them character- 
ized by three parameters: its redshift Zabs, its column density Nhi, 
and its Doppler parameter b = yjlkT jmn, where k is the Boltz- 
mann constant, T is the kinetic temperature of the gas and nin is 
the mass of the hydrogen atom. 

The redshift and column density characterising each absorber 
are drawn from a distribution of the form 

/(NHI,z) = yVo•(l+z)''■NH^^ (4) 

where Nq is a normalization constant. This function defines the 1- 
dimensional distribution of the Hi present in the IGM probed by a 
random line-of-sight. The number N abs of systems for each line-of- 
sight is drawn from a Poisson distribution with parameter 

{N,b.)= r r /(NHi,z)dNHidz, (5) 

where the integral is carried out over appropriate redshift- and col- 
umn density intervals /; and /nhh respectively. 
We use different sets of input distributions that include the evolu- 
tion of both low- and high density absorbers, and that give rise to 
the following models: 

MMC This model relies on the input distributions from IMadaul 
their equation 10) listed in Table [T] Here, the Doppler pa- 
rameter is kept constant at a value b = 35.0 km s"', which roughly 
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Tabl e 2. Typ es of absorbers and their corresponding parameters, adopted 
from lBCG99l their equation 10. 



Nhi [cm -2] 




No 


7 


/? 


10'2 - 10"* 




3.14- lO'' 


1.29 


1.46 


10'** - 1.59 


10 '7 


1.70- 10* 


3.10 


1.46 


1.59- 10" - 


1020 


1.90- 10** 


0.68 


1.50 



matches the mean fe- value derived by iRauch et alj ( Il99 2l). This 
model c orresponds to one of the several models presented by 
lBCG99l and is referred to as the MC-NH method in their work. 



BMC This model matches the best method of IBCG99| named MC- 
Kim. The corresponding parameters for the line-density evolution 
and column density distribution functions are summarized in Ta- 
ble |2] In this model, in contrast to the MMC model, the Doppler 
parameter for each absorber is drawn from a truncated, redshift- 
dependent Gaussian distribution of the form 



P{b) = &(b-b,,)- 



1 



V27! 



exp(-^(*-A')- 



where @{x) is the Heaviside function: 
®(") = 1 1, x>Q 

and the mean, standard deviation, and truncation value as a function 
of redshift are given by yu (z) = -3.85 z-h38 .9, o-(z) = -3.85 z+20.9 
and b„.{z) = -6.73 z -I- 39.5, respectively. [bCG^ originally used 
this model to analyse the impact of the intergalactic attenuation on 
high-redshift galaxy colors in the range 1.75 < z < 5.0, but we use 
it in the extended range 0.2 < z < 5.4 1. Since our highest redshift 
limit is not too far away from'BCG99''s, we may safely assume that 
the model is valid in our extended redshift range. 

MMC without Lyman limit systems In order to asses the impact 
of the Lyman limit systems on the intergalactic absorption, we 
introduce this model, which consists of the same input distributions 
as the MMC model, excluding the systems with column densities 
Nhi > 1-59 • 10 '^ cm -2. 

In all the models describe above, the attenuation factor for 
each absorber is given by exp[-T(/i)], where, for the general case, 
the absorption coefficient t {X) can be written as 



(6) 



The first term on the right-hand side is the opacity due to the ioni- 
sation of neutral hydrogen by photons with wavelengths /I ^ = 
91.18 nm. It is given by 

Tll« = Nhi •0az.z.-'i)-gW-o-»-|^j , (7) 

where iTo<, = 6.3 -10""* cm^ is the Hi photoionisation cross-section, 
and g is the Gaunt-factor for bound-free transition^ The second 
term is the sum of the opacities due to resonant scattering at each 



* An extensiv e tabulation of values for the Gaunt-factor can be found in 
iKarzaslfTgeil) . 



transition of the Lyman seriefl In general, the absorption coeffi- 
cient for the transition 1 i is 



Ti{X) = Nhi • o-j • 0(a,,x). 



(8) 



The cross-section cr, is a function of the Doppler parameter b, the 
oscillator strength of the transition ft, and the resonant wavelength 
/i,, and is given by 



(9) 



where A/lo = A, bic is the Doppler broadening, and the variable 
x = {X — /l,)/A/lo is the distance to the line center in Doppler units. 
We assume the profile function of the absorption line to be given 
by the Voigt-Hjerting function 



//(fl„x)=— jdy. 

J-oo (x - v)- + a] 



(10) 



Here, Uj = A^^Til{A7t/S.Ao) is the relative strength of the natu- 
ral broadening to Doppler broadening for the ;th transition, and 
y = w/b is the kinetic velocity in units of the Doppler parameter. In 
this work, we neglect the opacity due to the photoionisation term 
and consider only the first resonant transition, i.e. the Lyo- transi- 
tion, since this is the only one of interest in the wavelength range 
studied here. Furthermore, we use the approximation to H for val- 
ues of a and column densities characteristic for intergalactic Hi of 
lTepper-Garcla|j2006h . 

3.1 Tlie transmission factor O 

The cumulative absorption along a random line-of-sight of the flux 
fem of a source at redshift z«n is calculated according to expression 

fern ef») 



fobs (Aobs) ■ 



(11) 



1 +Ze 

where Aab.i and Ag„ are the observed and the emitted wavelengths, 
respectively. These are related by A„i,s = A^,,, • (1 + Zcm)- The quan- 
tity cD is the transmission factor and is given by 



O (^) = P] exp[-T2 (A/(l + z,)] = exp - J] T2 (A/(l + z,) 



,(12) 



where t is given by equation l[8]l, and z, is the redshift at the epoch 
of absorptior^. 

Introducing the relation fdA) = fg„, [A/{1 + z)]/(l + z), it 
follows from equations l|2). 1 11 U . and il2\ that 



1-£'a(z) 



— o 

J, 1, ■(!+,.) 



(A) dA . 



(13) 



The right-hand side of this expression is just wavelength-averaged 
value of O at redshift z, and we will denote it by O;. Note, how- 
ever, that this quantity still depends on redshift, as indicated by 
the subscript. Since Da (z) and (t>- differ only by a constant factor, 
they may be considered as equivalent with respect to their statistics, 
which will be discussed in the next section. 

2 We adopt the convention that the Lya transition (from the ground state to 
the next higher energy level) be identified with i = 2, the Ly/3 transition with 
i = 3, etc. The photoionisation cross-section is thus consistently denoted by 

O"o<j. 

^ The reader shall bear in mind that N hi , and <t , and a ,-through the depen- 
dence on the Doppler parameter-in equation (8) are dift'erent in general for 
each absorber. However, we do not write this explicitly by e.g. introducing 
a new subscript in order to avoid a cumbersome notation. 
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3.2 Distribution of Z) » 



Judging from the dependence of O on t and N ats (see equation 
112b . it is expected that the transmission factor and consequently Da 
are rather complicated random variables. Nevertheless, as is well 
known from statistics, a random variable x that can be expressed 
as the product of a large number of small, statistically independent 
factors is distributed lognormally, i.e. according to the distribution 



f(x;ij,cr) 



1 



(14) 



where fi and cr are the mean and the standard deviation of In x. This 
expression is equivalent to the statement that a random variable x 
is distributed lognormally, if and only if its logarithm is distributed 
normally. 

The expected value yu' and the standard deviation cr' of a log- 
normal distributed quantity can be expressed in terms of the param- 
eters fi and cr as 



(15) 



and 

/ 2 \l/2 

o-'=(e- -l) fj'. (16) 

From the form of equation ( 112b we may suspect that the transmis- 
sion factor is a lognormally distributed variable, since it can be ex- 
pressed as the product of a large number of statistically indepen- 
dent factors that take on values in the range [0, 1]. The implications 
of this statement are profound: If O is distributed lognormally, so 
does <I>; and consequently Da(z). Furthermore, due to the prop- 
erty of the lognormal distribution stated above, the effective optical 
depth o f the Lya absorpt ion, usually defined asr^ff = - \n(l - D a) 
(see e.ie. lKim et alj200lh should obey a Gaussian distribution. This 
result follows independently from the fact that the total optical 
depth can be expressed as the sum of the independent contribu- 
tion of each system, as indicated in equation ( 112b . Thus, for a suffi- 
ciently large number of absorbers Nab.,, and if the optical depth for 
each absorber has the same mean value (t) and dispersion cr (r) at 
each wavelength (i.e. redshift), Teff should obey a Gaussian distri- 
bution at each redshift, centered at A^^;,, (t) and with a dispersion 
■\/N^ cr (t). These statements are completely independent of the 
form of evolution of the intergalactic neutral hydrogen, as long as 
the transmission factor can be expressed in the form of equation 

(ED. 

In order to test whether the values of Da obey a lognormal 
distribution, we compute Da using the MMC and BMC models for 
an ensemble of A'los =4-10^ lines-of-sight at a given redshift, and 
from these values we compute the mean and standard deviation of 
In Da according to the equations 



A* (In Da): 



and 



1 



A'l. 



(T^dnDA) = V (IuD;, - (IuDa)) , 



(17) 



(18) 



With these parameters, we generate for each redshift a set of 
normally distributed random numbers using a standard random- 
number generator, and compare them statistically to the In Da val- 
ues from our simulations at each given redshift, in order to deter- 
mine whether the latter are normally distributed. Recall that In Da 
is distributed normally if and only if Da obeys a lognormal proba- 
bility distribution. As a first approach, we choose to compare both 



MMC (lognomial) 

MMC without LLSs (lognormal) 

MMC (Gauss) 

MMC without LLSs (Gauts) 




Figure 1. Values of d for the MMC model assuming a lognormal (heavy 
sohd line) and Gaussian (heavy dashed line) parent distribution. The shaded 
area indicates the range i5 ^ 0.05, that corresponds to a probabiUty of 1.9 ■ 
10"^ that the data are not drawn from the same parent distribution. Note how 
the assumption of a Gaussian pai'ent distribution improves with redshift, but 
that a lognormal distribution is a better assumption for low redshifts. The 
light solid and light dashed lines correspond to the MMC model without 
Lyman limit systems (see Section D.2. U . Note that the approximation to a 
lognormal distribution, is better when these systems are included (see text 
for discussion). The fluctuations seen with respect to a perfect smooth curve 
in each case are due to the random nature of the process at each redshift, 
and are not significant. 



data sets via a ;^'^-test similar method, and to this end we com- 
pute for each set of values the 50 per cent quantile (median), and 
the ±34.13, ±43.32, ±47.72, ±49.38, and ±49.87 per cent quantiles 
around the median. We define the following measure 



(19) 



where the 2's are the quantiles determined from the set of random 
numbers distributed normally, the ^'s are the corresponding quan- 
tiles of the In Da values at a given redshift that result from our 
simulations, and = 1 1 is the number of quantiles. By definition, 

5 explicitly gives the absolute deviation of one data set as charac- 
terised by these 1 1 quantiles with respect to the other data set, and 
thus quantifies the departure of the assumed distribution. Indeed, 
the smaller the value of S, the larger the probability that both data 
sets belong to the same parent distribution. The values of S com- 
puted in this way turn out to lie in the range (5 • 10"-, 5 • 10"'), 
which may seem vanishingly small. However, since the distribu- 
tion of the random variable 6 itself is unknown, we determine the 
significance of the results by generating two sets of normally dis- 
tributed random numbers with mean and standard deviations dis- 
tributed uniformly, and comparing them to the values found for our 
simulated data. We do this for 1.2 • 10* pairs of sets, and deter- 
mine from these what is the fraction of realisations with a value of 

6 smaller than a given value. According to these estimates, there is 
a probability of {1.4 • 10^^3.1 • 10-^ 1.9 • 10-^ 2.3 • lO^^l that 
a value 6 < {0.50, 0.25, 0.10, 0.05) will be produced by chance, 
respectively. In Figure [T] we show the values for 5 as a function 
of redshift for the full MMC model, assuming a lognormal parent 
distribution for Da (heavy solid line). The result are qualitatively 
the same for the BMC model and are not shown for clarity. Note 
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Figure 2. Distribution of tlie Z)a values at a given redshift computed from simulations based on the MMC model (histogram). The distribution is normalised 
to unit area. The bin size used to compute the histogram at each redshift has been arbitrarily chosen to be given by max[ Da (j)I/100, and is hence different 
at each redshift. The binning of the data, i.e. the Da values at each redshift obtained from the simulations, is intended only for display purposes. Shown are 
also the analytic lognormal distribution (solid line) computed according to equations )14t . )171 . and (18), and the analytic Gaussian distribution (dashed line) 
computed using the same equations as before but replacing InDA by Da- The y-axis indicates the probability of the corresponding Da value on the x-axis. 
Note that the binned data have not been used in any form to estimate the parameters of the corresponding probability distribution. Hence, the curves shown 
are not fits to the binned data, but are computed using only the mean and standard deviation of the raw data. Note the excellent agreement between the data 
and the analytic lognormal distribution, especially with respect to the skewness and the cut-off at D a = 0. The distribution of the data is not well described 
by a Gaussian distribution at the lowest redshifts shown. However, the approach of the data distribution to a Gaussian distribution increases with increasing 
redshift, as expected (see text for details). 



that for the whole redshift range shown, 6 < 0.5, and in particular 

6 < 0.1 for Zem > 1-5. 

It follows from the Central Limit Theorem that the larger the 
number of absorbers Nai,.i, the closer the approach of the distribu- 
tion of O to a lognormal distribution, since the distribution of In O 
approaches a Gaussian. Given that N„bs increases with redshift, it 
should be expected that the accuracy with which the distribution of 
<t> approaches a lognormal distribution also increases with redshift. 
Furthermore, if the integral in equation il3\ is approximated by a 
sum of the form 

(=1 

where is the number of pixels, it is apparent that for a sufficiently 
large A', and assuming that $ has the same mean value (O) and dis- 



persion it(0) at each pixel, the distribution of (S>- should approach a 
Gaussian distribution with mean (O) and dispersion cr(0)/A'. Since 
the width of a given restframe wavelength range, i.e. the number 
of pixels also increases with redshift as (1 -)- z), it is expected that 
the approximation of the distribution of O- and hence of Da to a 
Gaussian becomes better with increasing redshift. We thus have the 
superposition of two effects: On the one hand, the distribution of 
O at a fixed wavelength approaches a lognormal distribution at any 
given redshift, with the accuracy increasing with redshift. On the 
other hand, the distribution of O; and Da (z) approaches a Gaus- 
sian distribution with increasing redshift. The net result should be 
that (t>- and Da are distributed lognormally at low redshifts, and 
that their distribution approaches a Gaussian for higher redshifts. 
Since Da asymptotically converges to unity for very high redshifts, 
its distribution at these redshifts is expected to be highly peaked 
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around the mean value. This is naturally given by the fact that the 
dispersion of O- scales like cr(9)/N around (O), and that A' in- 
creases with redshift, as stated above. 

In order to test this expectation, we compute the analytical log- 
normal (Gaussian) probability distribution of Da as given by equa- 
tion l ll4b with the mean and standard deviation of IhDa (^^a) com- 
puted using equations l ll7t and l llSt from the ensemble of values 
that result from our simulations at each redshift. The comparison 
of the analytic distribution thus obtained to the distribution of sim- 
ulated Da values is shown for the full MMC model in Figure|2] The 
results for the BMC model are qualitatively the same and are hence 
not shown. Note the excellent agreement at all redshifts between the 
lognormal probability distribution and the distribution of the Da 
values resulting from our simulations. We want to emphasize that 
the curves shown are not fits to the binned data, but are computed 
using only the mean and standard deviation of the unbinned data. At 
lower redshift, the agreement at the lower cut-off, i.e. at Da = is 
worth mentioning. It is remarkable that this cut-off, which is phys- 
ically given by the fact that absorption as measured by Da can- 
not take on values smaller than zero, arises in a natural way due 
solely to the fact that Da is distributed lognormally. Note, in con- 
trast, that a Gaussian distribution does not satisfactorily describes 
the distribution of the data at these low redshifts; in particular, the 
Gaussian distribution does not display the sharp cut-ofl^ at Da = 0. 
Nevertheless, the description of the data by a Gaussian distribu- 
tion becomes better with increasing redshift, as expected, and also 
the corresponding Gaussian distribution becomes narrower at ev- 
ery increasing redshift (taking into account the change in the scale 
of the X-axis from the first to the last panel). Moreover, note that 
the analytic lognormal and Gaussian distribution become visually 
indistinguishable from each other at Zem > 3.0. In order to assess 
quantitatively the differences between these distributions with re- 
spect to the distribution of the data, we compute again the value of 
S at each redshift, assuming now that the Da values obtained from 
the MMC model are drawn from a Gaussian parent distribution. We 
compare these values to the corresponding values computed before 
for an assumed lognormal parent distribution. This comparison is 
shown in Figure [T] As can be seen, the values of d for both distri- 
butions are low at all redshifts, and the difference between them at 
Zem ^ 2 is negligible small. This explains why the lognormal and 
Gaussian distribution functions shown in Figure |2] are practically 
indistinguishable from each other. However, note that the values of 
S for a assumed Gaussian parent distribution become vanishingly 
small with increasing redshift, eventually becoming smaller than 
the corresponding values for the lognormal distribution. Besides, it 
can be see also that the 6 values for a assumed lognormal distribu- 
tion rise again towards high redshifts, implying that the distribution 
of the data at these redshifts is no longer well described by a log- 
normal distribution. This confirms our statement made above, that 
Da is expected to be distributed lognormally at low redshifts and 
normally at higher redshifts. The redshift at which the transition 
from a lognormal to a Gaussian distribution takes place may de- 
pend on the particular set of input distributions used. For redshifts 
where the distribution of Da is well approximated by a lognormal 
distribution, the optical depth should be distributed normal ly. This 
result is however not in agreement with t he res ults from e.g. lMadaul 
(I1995L his Figure I) iMeiksin & White! ( |2004 their Figure BI), or 
iBernardi et alj ( |2003| . their Appendix C), where the distribution of 
the mean flux, i.e. the averaged transmission factor O- is claimed 
to be highly inconsistent with a Gaussian distribution. It is usually 
assumed that the presence of the high-density Lyman limit systems 
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Figure 3. Curve of growth of the Lya absorption hne, for three typical val- 
ues of the Doppler parameter. The shaded region corresponds to the column 
densities characteristic to Lyman limit systems. 



is responsible for this behaviour, and we next make use of our sim- 
ulations to address this question. 

3.2.1 The Effect of the Lyman limit systems (LLSs) 

In order to quantify the effect that the LLSs have on the distribu- 
tion of Da, we compare, following the analysis of Section [T2] the 
results of the simulations for the MMC model with and without the 
optically thick LLSs. We find that Da is lognormally distributed as 
well with a high confidence for the case where the LLSs are ex- 
cluded, as can in Figure[T] where the values for 6 computed at each 
given redshift are shown. Again, the results for the BMC model 
are qualitatively the same and are omitted here. Note that the val- 
ues of 6 for the MMC model without LLSs are larger at redshifts 
Zem ^ 3 than those for the full MMC model. This is due to the 
fact that the absence of the optically thick LLSs enhances the prob- 
ability of the attenuation factor exp (-r) to be closer to unity at 
a given wavelength along a random line-of-sight. This means ef- 
fectively that the number of factors in equation JI2b are reduced. 
Hence, when the LLSs are absent, the approach to a lognormal dis- 
tribution should be worse with respect to the case where the LLSs 
are included. Furthermore, this effect should be enhanced towards 
lower redshifts, for which the number of factors, i.e. of absorbers 
decreases as oc (1 + ^)''. Note that the trend is the opposite in the 
case of a assumed Gaussian parent distribution. 

We hence find when computing the evolution of Da with 
the MMC model with and without LLSs that the predictions for 
the evolution of Da with redshift for the full MMC model are 
practically indistinguishable from the results of the MMC model 
without LLSs, which implies that the effect of these systems 
on the total absorption is small. Furthermore, it turns out that 
the distribution of Da is lognormal irrespective of the presence 
of LLSs. All these results point to the fact that the LLS have a 
negligible impact on the evolution of Da, as long as the input 
distributions used here correctly describe the number density 
evolution of the absorbers. This would confirm similar results 
already found bv Desiacqu es et al] j2007h . Also, as pointed out by 
jMcDonald et al. (2005), the picture would not be fundamentally 
different if higher column density systems such as damped Lya 
systems were present. We will come back to this point in Section 
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We think that, as long as the distribution functions realistically 
describe the evolution of the Lyo- absorbers, the LLSs cannot truly 
have a great impact neither on the absorption as measured by Da, 
nor on its statistics because of the following reasons: 1) They are 
scarce, and even more compared to the thinner Lyo- forest systems 
(e.g. of the order of 5 LLSs on average along a random line-of-sight 
compared to approximately 10^ Lyo- forest absorbers out to z = 3.0, 
according to the MMC model); 2) Their contribution to the absorp- 
tion is due solely to the absorption at the resonant wavelength of 
the Lyo- line, and the equivalent width of the Lya absorption line 
arising in systems with column densities Nhi = 10" -'""* cm"^ is 
not too different from that of the systems with column densities 
Nhi = lO'^^'^cm"- (cf. the curve-of-growth in Figure[3]l, which 
are by far more numerous. In order to estimate the relative contribu- 
tion of each population at a given redshift, we weight the equivalent 
width W(N, b) of the Lyo- absorption line as a function of the col- 
umn density N and Doppler parameter b with the column density 
distribution, and compute the ratio 



Q(Zei„\b) ; 



where 



hya 



(20) 



(W), = r™ r """ (1 + z) • Wo (N,b)f, (N,z)dNdz, (21) 

Jo J N,„i„ 

Here, /, is the distribution function of population / e {hya, LLS), 
Wo (N, b) is the rest equivalent width of a Lya absorption line for a 
column density A' and Doppler parameter b, and A'„„„ and N are 
the column density limits that define each population, respectively. 
For a reasonable value for the Doppler parameter of 36 km s "' , and 
the input distributions of the MMC model, we find that g » 1 at all 
redshifts. The result is qualitatively the same for the BMC model. 
Hence, It follows from these models that the Lytt forest systems 
dominate the absorption over the optically thicker Lyman limit sys- 
tems at all epochs. This would explain at least in part why the dif- 
ference between the predictions for Da from the MMC model with 
and without LLSs, is small. Also, it is consistent with the theoret- 
ical expectation that the distribution of Da should not be far from 
lognormal or Gaussian, with or without LLSs, since this only de- 
pends on the fact that the absorption factor be expressed in the form 
of equation l ll2t . and this is truly independent of the form of the in- 
put distributions, as stated previously. 



3.3 Scatter in Da 

We expect the intrinsic scatter in the absorption due to cosmic vari- 
ance to be strongest at lines-of-sight of middle length. At low red- 
shifts, both the number of thin Lya forest clouds and thick Lyman 
limit systems is small, and the addition of a few more does not 
change dramatically the amount of absorption. However, the mean 
number of Lya forest clouds and LLS increases with z, and thus 
the probability of encountering more or less systems than average 
increases as well. Correspondingly, the absorption increases and so 
does its scatter. At even higher redshifts, the number of Lyo- forest 
systems increases so dramatically and the absorption is so severe 
that the addition of more systems does not make any difi'erence 
neither to the absorption nor to the scatter. Thus, we should expect 
the stochastic effect, i.e. the scatter in absorption, to peak at some 
intermediate redshift z„„. 

We compute the intrinsic scatter at the cr level for the values 



MMC 

MMC (no LLS) 
BMC 



Figure 4. Evolution of the intrinsic scatter of D a due to the stochastic na- 
ture of the absorption in the intergalactic medium for the two competing 
modefs MMC (dashed hne) and BMC (dotted hne). The dot-dashed hne 
corresponds to the MMC without Lyman limit systems (see Section 13.2.11 
for details). Note that the behaviour of cr(Z)A) is qualitatively the same for 
all three models. As expected, the amount of scatter for the MMC model is 
less in the case where the Lyman limit systems are absent. 



of Da at each given redshift obtained from our simulations using 
equations l |15t , l |16t , l |17| l, and l |18l l. The result is shown in Figure 
m where we show the evolution of the scatter in Da with redshift 
for the MMC model with and without LLSs. We also include for 
comparison the result from the BMC model. It can be seen that, 
irrespective of the model, the scatter peaks at a intermediate redshift 
between Zem ~ 3.5 and Zem ~ 4.0. Note that the peak is significant, 
since it represents an increase in the scatter of 2.5 times with respect 
to its valu e a.iZem ~ 1 - It is interesting that this result had also been 
found by IZuol ( 1 19931 . their Figure 2), who using a semi-analytic 
approach and different input distributions, reported that (t(Da) is 
largest at redshifts near 3.7. Thus, the qualitative behaviour of the 
intrinsic scatter of Da shown in Figure |4] may be an unavoidable 
feature of this observable, which could explain at least in part the 
large scatter in the measurements of Da seen in the same redshift 
interval (cf. Figure [6ll. As will be shown in Section |5] however, 
there is a disagreement at the 3 a level between the amplitude of 
the scatter in the observations at these redshifts and that predicted 
by the models. 

When comparing models that only differ by the presence of 
the optically thick Lyman limit systems, we find that the scatter in 
Da is larger at any given redshift when the LLS are present, as 
expected. None the less, the absolute value of the scatter does not 
differ significantly between the situation where these systems are 
present and where they are absent. Thus, the net effect of the LLSs 
is to enlarge the intrinsic scatter in the absorption, contributing only 
marginally to the mean value of the absorption itself. However, note 
that for redshifts Zem ^ 4.0, the amount of scatter at the Icr level 
becomes indistinguishable between both scenarios, owing to the 
fact that the thinner Lyo- forest systems are overwhelmingly more 
numerous at these redshifts. 

The results stated above are consistent with the previous con- 
clusion obtained from our simulations that the optically thick LLSs 
cannot have a great impact on the absorption as measured by Da 
and on its scatter. 
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4 OUR Oa measurements 

We want to compare the predictions for tlie evolution of £)a 
tliat result from the models described above to observations. For 
this purpose, we use previous measurements of Da reported in 
the literature, and we perform ourselves new measurements of 
this quantity, using QSO spectra from the SDSS Data Release 5 
([Adelman-McCarthv et al. 2007), thus extending the redshift range 
of the measurements to Zem = 5.41. 

Our selection procedure of sources suitable for this purpose 
was as follows: Since the wavelength range available from SDSS 
(DR5) isAe [380, 920] nm, and we measure the continuum depres- 
sion £)a in the restframe wavelength interval /I 6 [105, 117] nm, the 
redshift of our sample is restricted to z,„,„ > 380/105 - 1 = 2.62. 
We choose z„,;„ = 2.7 as our lowest redshift in order to avoid the 
low S/N at the blue end of the spectrograph. A simple query for 
high-z quasi-stellar objects on the SDSS SkyServer Spectroscopic 
Query Form with this restriction alone returns around 2400 spectra. 
From this first selection, we rejected those objects for which the 
redshift was either not measured, the measurement had failed, or 
the measured photometric and spectroscopic redshifts were incon- 
sistent with each other. We binned the quasars in redshift intervals 
of Az = 0. 1 and selected for each redshift bin the spectrum with the 
highest S/N, leaving us with 28 sources, from which we removed 
three further objects due to low data quality. The resulting sample 
is listed in Table[3] 



4.1 Continuum fit 

The continuum of a quasa r is often assumed to be of the 
form /, = • v"' (see e.g. ISteidellll987bl , iLaor et alj[T997h . or 
equivalently, /,i = /]" • /l""'', where both indices are related by 
= 2- ay, and /]• = fy • c'""", and c is the speed of light. Em- 
pirical evidence in favour of the assumption of an u nderlying con- 
tinuum in form of a power-law has been given by e.g. lO'Brien et al.l 
( Il988 ^. who derived an average power-law index of o-^ = 2.36 for 
A 6 [80, 121.6) nm and ay = 0.67 for A e [12 1.6, 190] nm. This 
is con sistent with the result already found by iNeugebauer et al.l 
that the spectral index a^ varies over large wavelength 
ranges. Izheng et al] ( 1 19971) found by constructing a composite QSO 
spectrum from 284 HST FOS spectra that a single power-law de- 
scribes well the continuum for wavelengths between 105 and 220 
nm with a V = 0.99 ± 0.05, but that the continuum steepens (flattens 
in A) significantly for A < 105 nm (a,, = 1.96 ± 0.15). Similarly, 
iTelfer et al. (2002) reported, using a sample nearly twice as large 
as the lZhe ng et al.'s sample, that the continuum in the extreme UV 
region between A = 50 nm and A = 120nm and in the near UV 
between A = 120nm and A = 300 nm is well described by a sin- 
gle power-law with Uy = 1.76 ± 0.12 and d,, = 0.69 ± 0.06, re- 
spectively. Furthermore, by constructing a composite QSO spec- 
trum from a homogeneo us sample of over 2200 SDSS QSOs, 
IVanden Berk et alj ( 1200 ll) concluded that the continuum in the rest- 
frame wavelength range 121.6 < A < 500 nm can be very well 
modeled by a single power-law either in wavelength or frequency, 
with Uy = 0.44 ±0.1. However, they did not estimate a spectral 
index for wavelengths A < 121.6 nm. While it is the consensus that 
the optica l/NIR continuum of quasar spectra is due to synchrotron 
emission jWhiting et al.|[200lh . the physical origin of the quasar 
UV/optical continuum is not yet established. On the hand hand, the 
observed soft X-ray/UV/optical spectral shape of quasars is found 
to be consistent with fre e-free emission o f optically thin gas at tem- 
peratures 10' - 10*^ K tearvainis|[l993h . which is well described 



by a power-law. On the other hand, it has been argued that the 
emission from a geometrically thin, optically thick accretion disk 
also reproduces well the obs erved quasar SED in this spectral re- 
gion dKawaguchi et al.|[200lh . Either model, however, is consistent 
with an underlying continuum in the form of a power-law. Alter- 
natively, as already mentioned in Section[2l other authors prefer to 
fit quasar continua locally using spline functio ns (see e.g. ILu et al.l 
1 19961 iKirkman et aljliool iTvtler et alj|2004ij|bl) searching for re- 
gions apparently free of absorption blueward of the hya emission 
line. 

Since the Lyor forest region is severely absorbed due to inter- 
vening Hi systems, especially for high-z QSOs, and a local fit to the 
continuum in this regi on is difficult, with the uncerta inty increasing 
with redshift (see e.g. iFaucher-Giguere et al.ll2007h . we choose to 
estimate the continuum of our selected sources in the hya forest 
region by fitting a power-law to the QSO spectrum redward of the 
Lyo- emission line and extrapolating it for A < 121.567 nm. As long 
as the assumption of the underlying power-law holds, this approach 
has the advantage that the continuum estimate is completely inde- 
pendent of the spectral resolution and S/N in the hya forest region. 
Nevertheless, this method may lead to an systematic overestimate 
of the continuum flux and thus to an corresponding overestimate 
of Da. It follows from the results quoted above that the quasar's 
intrinsic continuum in spectral regions below and above the Lyo- 
emission line may both be well described by power-laws, but with 
potentially different spectral indices. More specifically, the spec- 
trum in the region below 121.6 nm (rest-frame) may be steeper in 
frequency, and hence flatter in wavelengtlfl with respect to the con- 
tinuum in the region above 121.6 nm (see e.g. Telfer et al, 20o3. 
their Figure 4). Thus, fitting a continuum in the form of a power- 
law at wavelengths A > 121.6 nm and extrapolating it towards 
A < 121.6 nm may well lead to an overesti mate of the continuum 
level in that region, as first pointed out by Seljak et alj j2003l) . ac- 
cording to whom the Da may be overestimated by at least 0.05 at 
Zem = 2.72. We will take this systematic bias into account when es- 
timating the uncertainty in our Da measurements (cf. Section |42l (. 

The power-law that we fit to each spectrum is of the form 

fAA;a,) = f°-(A + Ao)-'". (22) 

where the flux amplitude the wavelength off-set Aq, and the 
spectral index a^ are the parameters to be determined. We fit the 
continuum in the wavelength range [130 ■ (1 + Zem), 900] nm, in or- 
der to avoid the red emission wing of the Lya line and the red end of 
the spectrograph, respectively. Note that the available wavelength 
range decreases with redshift, and this may introduce an increasing 
uncertainty in the fitted continuum. The fit parameters and their un- 
certainties were obtained with help of the IDL task CURVEFIT. It 
turns out that the continuum fit is rather insensitive to the uncertain- 
ties in the flux amplitude and the wavelength off-set, and extremely 
sensitive to the uncertainty in the spectral index. Because of this 
and for simplicity, in our further analysis we neglect the error in 
the first two parameters and consider only the uncertainty in the 
spectral index to be of relevanc^fl An example of a QSO spectrum 
and its corresponding fit are shown in Figure[5] The spectral index 
and its uncertainty for each source are listed in Table [3] (columns 
three and four, respectively). The uncertainty in the spectral index 



Recall that o;^ = 2 - Uy, which implies da^ = -day. 
^ The full list of fit parameters for each source and 
their uncertainty are available in machine-readable form at 
|http://astro.physik.uni-goettingen.de/~tepper/da/fitparam.txt| 
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Table 3. QSO sample selected from the SDSS DR5. The first four columns include the object designation", the emission redshift as quoted in the SDSS DR5 
catalog, the spectral index used for fitting the continuum, and its uncertainty, respectively. The last three columns list the measurement of D a and its total 
uncertainty, and the uncertainty due to different error sources, respectively (see text for details.) 

Object Zem a A O-fl, (ff,!) Da O-flt(DA) crsys(£>A) 
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" The designation of each object meets the lAU nomenclature, as requked. 
www.sdss.org/dr5/coverage/IAU.html 

quoted in column four and denoted by tTd, (o- ,) includes the for- 
mal error inherent to the fitting process, which is due to the error in 
the detected flux, the removal of broad emission lines in the region 
/i > 121.6 nm - which makes the wavelength range become patchy 
and the systematic uncertainty due to the decreasing available 
wavelength range with redshift mentioned above. The systematic 
efi"ect of the latter in particular can clearly be seen on the overall 
trend for the uncertainty in the spectral index to increase with red- 
shift. 
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For details on the official SDSS designation of an object, please consult 



4.2 Measurement of D a 

For each QSO spectrum, we compute the total absorption pixel by 
pixel in the restframe range A e [105, 117]nm according to the 
expressior^ 

£'A(z™;Q'i) = 1 - -T— > , -77- -, (23) 

where A'p,> is the total number of pixels between 
/i] = 105 • (1 + ZOT,)nm and A2 = i\l ■ (\ + Ze,„)nm. The lower 
limit of this range is chosen as to avoid the wings of the Lyy6 + 
Ovi emission lines, while the upper limits is such that the blue 
wing of the Lya emission line is excluded in the computation of 



* This is just the discrete version of equation (2). 
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X [nm] 

Figure 5. Spectrum of the QSO SDSS Jl 12107.99+5 13005.4 at 
Zem = 2.843. The heavy solid and heavy dashed line indicate the power- 
law fit to the continuum and its uncertainty, respectively. The corresponding 
spectral index is = 1.256 ± 0.003. 



Da- The measurement of Da for each source together with its 
corresponding uncertainty are hsted in Table |3] (column five). 

There are two independent sources of uncertainty in the mea- 
surement of Da^ 1) The error in the spectral index inherent to the 
continuum fitting, i.e. cto, (a^), and 2) the possibility that the con- 
tinuum level may have been systematically overestimated due to the 
change in slope of the power-law atA= 121.6 nm (cf. Section|4T|l- 
We shall refer to the uncertainty in Da due to the first and second 
sources crat (Da) and (TsysiDA), respectively. The total uncertainty 
in the measurement of Da -quoted in column five of Table [Sj- is 
then assumed to be given by cr^ (Da) = c^^^ (Da) + cr^ys (Da)- 

As a consistency check, we adopt two different 
methods to estimate (Tf,t(DA). The first method assumes 
that the ±cr range for each measurement is given by 
-Da [z«„; ff,! ± o"!,! (a: a)] -Da U™; tt^], respectively, for the 
corresponding values of o"fit (a a) listed in Table [3] The second 
method is based on error propagation, according to which the error 
cr(/) in the estimate of a quantity /(x,), which depends on n 
independent random variables {x, j, each with an uncertainty cr(x,), 
is given by 



the values computed according to method 1 , which are listed in Ta- 
ble[3](column six). 

It is not clear based on the results found in the literature (cf. 
Section |4jTJ whether the continuum blueward of Lyo- may always 
be overestimated when assuming it to be described by an extrapo- 
lation of the power-law redward of hya. However, we will assume 
that this is the case in order to allow for the possibility that our 
Da measurements may have been systematically overestimated. In 
order to compute the corresponding uncertainty in Da, we make 
the following assumptions: First, the continuum in the region blue- 
ward of Lyff shall be described by a power-law as well, but with a 
smaller spectral index - i.e. the continuum is flatter - with respect 
to the region redward of Lya. Since it is not possible to fit a power- 
law in this regime for the data we use, we assume that both power- 
laws (above and below Lyo-) are normalised at A = 121.6 nm. Fur- 
thermore, and for simplicity, we take the difference in the spectral 
indices below and above Lyo- to be equal 1 (which corresponds 
roughly to the findings in the literature quoted in Section l4m i. Un- 
der these assumptions, the value of Da for our sample turns out to 
be overestimated on average by less than ten per cent. This value 
is slightly lower for our measurement of Da at Zem = 2.71 than 
the correction estimated bv lSeliaketai] ( l2003h . However, they also 
remark that the exact value of the bias depends on the particular 
method applied to extrapolate the continuum in the Lyo- Forest re- 
gion. We hence make the rather conservative assumption that each 
measurement of Da has an additional uncertainty cTsys (Da) of ten 
per cent, and add this uncertainty only to the low error bound, since 
we are assuming that the continuum, and hence Da, may have been 
overestimated. We could in principle correct our Da estimates by 
shifting all measurements by ten per cent towards a lower value, 
but since it is not clear whether the continuum has been overesti- 
mated for every single measurement, we prefer to express this bias 
as an uncertainty. The uncertainty cr^ys (Da) for each measurement 
is listed in the last column of Table|3] 

Note that we do not correct our measurements for contamina- 
tion of metal lines. However, this should not in troduce a large error, 
since their contribution is small. For example. iTvtler e t al. (2004^ 
find that they contribute by 2.3 ± 0.5 per cent to the total absorption 
at z = 1.9. The validity of this assumption will be tested in the next 
section by the comparison of Da -measurements to the results of 
our simulations, in which the absorption due to metal lines is not 
included. 



o-^(/) = |;.V,).^ 



(24) 



In the case of our measurements of Da, only the uncertainty in the 
index a a is relevant, and hence equation ( 124b becomes, with the 
appropriate notation, 



,(Da) = o-fi,((!rj)- — y ln(.l, + ^o) 



/ofa(^i) 

fc(Ai;aA)' 



(25) 



where we have used equations i22i and ( 1231 ). Note that we do not 
include the term due to the error in the detected flux fohs, since this 
is already included in crot (a a). The error computed according to 
the first method is asymmetric, as expected, while the error com- 
puted using equation J25b approximately (in some cases exactly) 
corresponds to the arithmetic mean of the former. Since the error 
estimates computed according to both these methods are consistent 
with each other, we may use either to estimate the uncertainty in 
Da due to the uncertainty in the continuum fit; we choose to use 



5 RESULTS & DISCUSSION 
5.1 Observations vs. Models 

We compute the evolution of Da in the redshift interval 0.35 < 
Zem < 6.0 using the models MMC (with and without Lyman limit 
systems) and BMC presented in the Section[3] which include the ef- 
fect of different populations of absorbers, i.e. Lya forest clouds and 
Lyman limit absorbers. For comparison , we include also the values 
of Da computed using iMeiksinI ( l2006 l)'s model (see Section [TTTT l. 
We refer to this model as MTC. For the MMC (with and without 
LLS) and BMC models, we simulate an ensemble of A'los = 4- 10^ 
lines-of-sight at fixed redshift, and compute for each of them the 
flux decrement Da according to equation i23l . In this way we get 
for each given redshift an ensemble of an equal number of Da val- 
ues for each model, from which we estimate the 50 per cent quan- 
tile (median), and the ±34. 13, ±43.32, ±47.72, ±49.38, and ±49.87 
per cent quantiles around the median, which correspond to the 
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±1, ±1.5, ±2, ±2.5 and ±3 cr ranges. We do not compute the mean 
and cr ranges in the standard way, since the distribution of Da is 
strictly speaking unknown a priori. However, according to the re- 
sults from Section [3] the distribution of Da is not too far from a 
lognormal or even a Gaussian distribution; hence, the identifica- 
tion of mean with median, of ±cr range with the ±34.13 quantile 
around the median, and so on, is justified. Since the MTC model 
only provides a (mean) transmission function at each redshift, we 
compute for this model one single value for Da at each redshift by 
numerically integrating the corresponding transmission function in 
the restframe wavelength range [105, 117] nm. 

We compare our simulations to measurements of Da accu- 
mulated in the literature over the past two decades, performed with 
different methods and approaches. This compilation is by no means 
intended to be complete. The reason for choosing these measure- 
ments is mainly that they were performed in more or less mutu- 
ally exclusive redshift ranges which all together cover the range 
< Zcm < 4. Our measurements extend this redshift range out to 
Zem = 5.41. Therefore, these previous measurements together with 
our own allow a comparison to the models presented here over a 
wide redshift range, but with the caveat that using measurements 
from difi^erent groups may introduce an artificial bias in the evo- 
lution of Da due to the differences in the particular method used 
and the different redshift ranges probed in each case, as previously 
mentioned. The literature data and our measurements are shown to- 
gether with our model calculations in Figure|6](see cited references 
for details on the corresponding measurements). 

As can be seen in this figure, the predictions for the evolution 
of the mean Da from all three models, MMC, BMC, and MTC, are 
practically indistinguishable from each other when compared to ob- 
servations for redshifts z„„ < 3. We do not include in this figure the 
predictions for Da based on the MMC model without LLSs, since 
the difference between this and the full MMC model is negligi- 
ble. At higher redshifts, the values of Da obtained from the MMC 
and MTC models match the observations quite well. Note that in 
spite of the fundamentally different approaches and input physics 
of these models, their results are very similar. As a result, it is not 
possible yet to discriminate between these models due to the uncer- 
tainty and the strong scatter in the observations, especially around 
Zem ~ 3.5 (cf. Section[33}. In contrast to the MMC and the MTC 
models, the predictions from the BMC model strongly deviate from 
the measurements of Da at z„„ > 3, underestimating the absorp- 
tion. Since the MMC and BMC models differ only in terms of the 
input distributions and not in terms of the method, this result sug- 
gests that the evolution of the Lyo- absorbers implied by the input 
distributions of the BMC model is slower than expected from the 
evolution of Da (cf. equation[3]l. If used to compute the magnitude 
changes for hig h-redshif t galaxies due to intergalactic absorption 
as done by e.g. lBCG99l the difference in the predicted evolution 
of Da between the MMC and BMC models corresponds to a dif- 
ference of slightly more than 0.6 mag in the predicted magnitude 
change for e.g. an Sd-type galaxy -with a constant Star Formation 
Rate- at z„„ = 4.0 in the F450W filter. The corresponding differ- 
ence between the MMC and MTC model amounts to less than 0.3 
mag. 

It is wo rth mentio ning at this point that, in contrast to previous 
models (e.g. lZuoll993h , we explicitly avoid normalising in any way 
the used distributions to match the observed Da at some given red- 
shift, or manipulating the models whatsoever to satisfy any other 
restriction. We simply take the distributions as reported in the liter- 
ature, where they were determined directly by means of line statis- 
tics by the authors and references therein. We mention this in order 




z 



Figure 6. Evolution of the mean Z)a computed according to the MMC 
(dashed line), BMC (dotted line), and MTC (dot-dashed line) models, com- 
pared to observations performed o ver the pa s t twen ty ye ars by different 
groups and different methods: (1) Zuo & Lu h99^. (2) ISchneider et alj 
1 1991), (3) Kirkman et al. (2007), (4) O'Brien etalj jl988h . The data points 
display lo" error bai's. Note that, in spite of the heterogeneity of the ap- 
proaches to measure Z)a, its mean evolution as computed using the MMC 
and MTC models matches well the observations over the entire redshift 
range shown, while they disagree strongly with the BMC model at redshifts 
Zim > 3.0. Below this redshift, the models are practically indistinguishable 
from each other For c ompleteness, we inc lude the empirical fits to the ev o- 
lution of Da from: (5) lzhang etd]il997h . and (6) lKirkman et alj j2005l) . 



to emphasize the rather good agreement between the observations 
and the evolution of the mean Da predicted by all three models at 
Zem < 3, and by the MMC model at higher redshifts, further re- 
inforced by the similarity between the results of the MTC and the 
MMC models, which are completely independent from each other. 

We include in Figure|6]the empirical fits of Zhang et al. ( 199^ 
of the form Da (z) = D° e "<'+'' with D» = 0.01 and a = 0.75 
and of iKirkman eTai] ( |2005|) of the form Da(z) = A{1 + zV with 
A = 0.0062 and y = 2.75. These empirical fits match the observa- 
tions quite well at z„„ < 4.5, and may be useful as a rough esti- 
mate of the mean absorption at those redshifts. However, as stated 
before, they should be taken with caution, especially at high red- 
shifts. While Da asymptotically converges to unity as the redshift 
approaches the epoch of reionisation, z„,o„, these empirical fits di- 
verge. In contrast, the predicted evolution of Da from all models 
described above does satisfy the expected and observed asymptotic 
behaviour. The rate of convergence to this limiting value certainly 
depends on the particular set of input distributions used, as shown 
by the comparison between the MMC to the BMC model. In gen- 
eral, the stronger the number density evolution, the faster the con- 
vergence. In other words, different input distribution functions im- 
ply different values for z„,o„, and this information may be used as 
a further constraint on the accuracy of a particular set of input dis- 
tributions. Extending the computations of the evolution of Da with 
the MMC model to redshifts z„„ > 6.0, it turns out that Da is al- 
most unity at Ze m ~ 7.0. This valu e is slightly higher than the value 
of 6.5 quoted bv lFan et al.l ( |2002|) for the epoch of reionisation. 

The compilation of Da measurements shown in Figure |6] 
displays a large scatter at all redshifts, in particular at redshifts 
Zem > 3. Previously, concerns were raised about the possibility that 
the amplitude of the intrinsic scatter due to cosmic variance may 
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Figure 7. Evolution of the mean D a according to the MMC model (solid 
line) in the range 0.2 < Zem < (>■ Here we show again our measurements 
of Da and those fr omlO'B rien et a l] il988l) , together with the most recent 
measurements from lKijkmaieralJ j2007l) ^he data points display Icr eiTor 
bars. The white and shaded areas around the solid Hne indicate the intrin- 
sic scatter at the ± 1 and ±3 o" level, respectively, due to variation in the 
absorption from one random line-of-sight to another. 



be enhanced by combining measurements performed with different 
methods and using heterogeneous data samples. Hence, we select a 
subset from our compilation following the criteria stated below, and 
compare this observations to the scatter as predicted by our simu- 
lations. The results of this comparison are summarised in Figure|7] 
Here, we show the redshift evolution of Da as computed from the 
MMC model. We include the 50 per cent quantile, i.e. the median, 
and the ±34. 13, and ±49.87 per cent quantiles, which correspond to 
the ±1-, and ±3 cr ranges around the median at each redshift. Again, 
the difference between the predictions of the MMC with and with- 
out LLSs is negligible, and thus we show only the results for the 
full MMC model. We do not include the BMC model in this anal- 
ysis in view of the disagreement with the observations at redshifts 
Zem > 3. Our new data subset con sists of our own mea surements 
together with the measurements of lO'Brien et alj ([l_988'), who use 
power-law fits to estima te the continuum level of the quasars, and of 
the measurements from lKirkman et ^ j2007h , who fit the contin- 
uum locally using b-splines. We s elect these measurements for the 
following reasons: lO'Brien et al.l 's data allow a more direct com- 
parison to our data due to the similarity in method to estimate D^; 
iKirkmsm et alj l2007) 's data are the most recent and accurate mea- 
surements of Da available to date. Furthermore, these data sets in- 
dividually cover the largest redshift ranges, and together they con- 
stitute the smallest set of individual measurements fully spanning 
the redshift range 0.2 < Ze,„ < 5.41 with some overlap. 

Within the quoted uncertainty, the observations are well 
matched by the evolution of Da as predicted by the MMC model 
in the redshift range 0.2 < Ze,,, < 5.41, with the exception of some 
outliers around z„„ « 3, where the scatter of our observations 
alone as well as the scatter in the full sample of measurements 
is large, as can be seen in Figure |6] It is apparent from Figure |7] 
that the largest scatter in the observed values of Da at a given red- 
shift cannot be ascribed to the intrinsic scatter in the absorption 
due to variations from one line-of-sight to another as computed 
with the MMC model, not even at the ±3cr level. Ne vertheless, 
note that this is only the case for lO'Brien et alj jl988h 's and our 



measurements, while lKirkman et al] ( l2007l) 's measurements are all 
nicely contained within the ±3 <r envelope of the model. This sug- 
gests that the observed scatter -or at least part of it- is not real, 
but only an artifact introduced by combining observations based 
on difi'erent methods to measure Da. This assumption is supported 
by the small amplitude of the scatter in the dist ribution of high- 
redshift galaxy colors already found bv IbCG9^ (cf. their Figure 
6) for a variety of inp ut distribution functions. Also, the fact that 
iKirkman et alj bOOTP 's measurements are well matched by the 
M MC model points to the fact that the particular method chosen 
by lO'Brien et alj ( 1 1988) and ourselves to estimate the continuum 
in the hya Forest region does introduce a bias in the measurements 
of Da. On the one hand, there is very likely a bias do to the al- 
ready discussed syst ematic overestimate o f the continuum level. 
Furthermore, neither [o'Brien et" d] jl988t) nor we corrected our 
corresponding measurements for metal-line absorption. These facts 
together would explain in part why precisely our and their measure- 
ments are higher on average than the median value of Da predicted 
by the MMC model, as can be seen in Figure[7] Even though its ef- 
fect is expected to be small, including the absorption of metal lines 
should eventually increase the agreement between the results from 
the models and the observations. 

There is thus rather strong evidence that most of the ampli- 
tude in the observed scatter in Da displayed in Figures |6] and |7] 
is introduced by combining measurements performed on heteroge- 
neous data samples using different methods. If real, however, the 
observed scatter would indicate that the models, in particular the 
MMC model, cannot account for the variation of Da among difi'er- 
ent lines-of-sight, even though it reproduces well its mean redshift 
evolution. 

According to the results from our simulations (cf. Section[3]2j, 
it is expected that Da as predicted by the MMC model be lognor- 
mally distributed with a high confidence at redshifts Zem i 6, and 
consequently, that In Da and hence the optical depth of Hi should 
obey a Gaussian distribution. However, this result is not supported 
by the observations given the poor agreement between the scat- 
ter predicted by the MMC models and the observed distribution 
of Da at a given redshift. On the other hand, the expectation that 
Da should be either lognormally or Gaussian distributed is based 
solely on the fact that the transmission factor be written in the form 
of equation l ll2t . which is certainly independent of the particular 
model accounting for the evolution of the hya absorbers. A dif- 
ferent set of input distribution functions may have an efi"ect on the 
mean value and the spread of the distribution of Da at each redshift, 
but not the on the form of the distribution. It may be thus possible 
that a difi'erent (more accurate) set of input distribution functions of 
the absorber properties may help to reconcile the lack of agreement 
between the amplitude of the predicted and the observed scatter, 
but even in this case the disagreement regarding th e form of the 
distribution when compared to previous s tudies {e.g. 



Madaij ll99 4 
White! l2004i) 



iBernardi et al J l200l iTvtler et alj l2004al iMeiksin & White! 
remains. 

Even though the amplitude of the observed and simulated scat- 
ter do not quite match, it is interesting to note that the redshift 
at which the intrinsic scatter in evolution of Da peaks (cf. Fig- 
ure |4ll roughly coincides with the redshift at which the observa- 
tions show a stro ng scatter (cf. Figure[6) . This redshift matches the 
redshift at which ISernardi et aU ( l2003h reported a local decrement 
in the evolution of the Ly o- optical depth Teff, recently confirmed 
by Faucher-Giguere et alj l l2007i) . which has been interpreted by 
iBemardi et alj as a signature of the reionisation of He ii. Note, how- 
ever, that at this redshift the Da measurements lie predominantly 
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above the mean value as given by the models (cf. Figure (6)1, and 
hence the trend suggested by these observations is exactly the op- 
posite from that found bv lBemardi et alj l l2003ll . since a decrease 
in Teff implies a decrease in Da- This strongly indicates that the 
difference in method and the heterogeneity of the samples of the 
different surveys does introduce an artificial bias in the observed 
evolution of £)a, in particular in its scatter. 

We finally want to highlight the following curiosity: By taking 
a close look at Figures |4] and [7] it becomes apparent that the max- 
imum in the evolution of (t(£)a) roughly coincides with the point 
of inflection of the curve that describes the expected redshift evo- 
lution of the mean Da- Mathematically, this would imply that the 
intrinsic scatter of £)a is proportional to the rate of change of the 
mean Da with redshift, i.e. 

<r(DA)oc— . (26) 

oz 

Indeed, if one computes numerically the derivative of {D^} with 
respect to z, it turns out that it qualitatively matches the evolution 
of (t(Da), up to a scale transformation. A physical interpretation 
may be gained in light of equation l|3): A rapid evolution of the ab- 
sorbers, quantified by the parameter y, implies a stronger evolution 
of Da with z. Since the evolution is however different in general 
for different lines-of-sight, this in turn implies a larger variation in 
the absorption from line-of-sight to line-of-sight, and thus a larger 
value of a (Da). This is mathematically consistent with the fact 
that the scatter given by the derivative of equation Q with respect 
to z is proportional to y. 



6 SUMMARY & CONCLUSIONS 

(i) We measured the cosmic flux decrement Da and its uncer- 
tainty due to statistical and systematic errors for 25 QSOs of the 
SDSS DR5 catalog in the redshift range 2.71 < z„„ < 5.41. 

(ii) We modeled the redshift evolution of the mean Da and its 
distribution at each given redshift in a Monte Carlo fashion, adopt- 
ing tw o among the various models presented by iBershadv et al.l 
( Il999h . We found that the predictions of the MMC model for the 
evolution of the mean Da reproduce well the observations in the 
range 0.2 < Zem < 5.41, in contrast to the BMC model. We con- 
clude from this that the underlying input distributions of the BMC 
model may not accurately account for the evolution of the Lyo- ab- 
sorption on the IGM. Hence, estimates of the impact of the inter- 
galactic attenuation on the photometric properties of high-redshift 
galaxies using this particular model should be taken with caution. 
Incidentally, by showing the rather good agreement between one 
particular set of simulations and the data, we show the power of the 
relatively simple approach used here to model the effect of inter- 
galactic absorption, as compared to models based on hydrodynam- 
ical simulations which are by far more complex. 

(iii) Through the comparison between the MMC and BMC mod- 
els, we showed that dilferent input distribution functions may have 
very different predictions regarding the evolution of Da, a quan- 
tity which is intimately related -and hence very sensitive- to the 
amount of Hi present along the line-of-sight. Therefore, we state 
that any model of the intergalactic absorption should first of all re- 
produce the observed evolution of Da and its distribution before it 
is used to e.g. correct synthetic or observed spectra for intergalactic 
absorption. 

(iv) The results from our simulations suggest that the distribu- 
tion of Da at a given redshift be well described by a lognormal 



distribution at low redshifts and by a Gaussian distribution at high 
redshifts, in agreement with our theoretical expectation based on 
the fact that the absorption is mathematically expressed as the prod- 
uct of small, statistically independent factors. This result implies 
that at redshifts where Da is distributed lognormally, the effective 
optical depth of the intergalactic Hi should obey a normal distribu- 
tion, contrary to the results of previous studies. However, this result 
should be taken with caution in light of the fact that the models, in 
particular the MMC model, cannot reproduce the amplitude of the 
observed scatter in Da. 

(v) We argue that most of the observed scatter in Da is intro- 
duced by combining measurements based on different methods. 
Nevertheless, the reason for the lack of agreement between the 
scatter predicted by the MMC model and the observed distribu- 
tion of Da is not ultimately settled. Because of this and until this 
discrepancy is clarified, we may warn about using any attenuation 
model based on the input distributions of the MMC model to esti- 
mate magnitude changes in the spectra of background sources. 

(vi) A larger, homogeneous sample of accurate measurements of 
Da over a wide redshift range is needed in order to allow for a more 
faithful comparison to models, and in particular, to determine the 
intrinsic /orm and amplitude of the distribution of Da as a function 
of redshift. 
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